Elastic solution of surface loaded layer with couple and surface stress effects

In this study, an elastic solution of an axisymmetrically surface-loaded thin layer resting on a rigid substrate is established by taking the surface stress and material microstructural effects into account. Derived solutions provide not only a means to investigate the size effects on the mechanical response but also a set of fundamental solutions essential for tackling contact problems in a micro/nano scale. In the formulation, the couple stress and surface elasticity theories are adopted to simulate the microstructured bulk layer and the surface material, respectively. A general solution of an elastic field within the bulk layer is obtained first by Hankel transform method and subsequently used together with the surface equations and boundary conditions to form a set of conditions essential for determining all unknown constants. After being fully tested with available benchmark solutions, results are used to study the role of surface and couple stresses on the load transferring mechanism to the substrate and its size-dependent characteristic for a wide range of external length scales relative to the internal length scales.

www.nature.com/scientificreports/ region of radius a and free of traction elsewhere. In the formulation and solution scheme presented further below, a reference cylindrical coordinate system {O; r, θ, z} with the origin O located at the center of the loading region, the r-axis directing along the infinite direction of the layer, and the z-axis directing downward is employed.

Field equations for bulk part.
To simulate the elastic response of the bulk material with microstructures, a fundamental couple stress theory proposed by Mindlin and Tiersten 19 and Koiter 23 is adopted. Basic equations (i.e., equilibrium equations, constitutive laws, and kinematics) governing the elastic field under the axisymmetric deformation and zero body force and couple are given by 66,67 where {σ rr , σ θθ , σ zz , σ rz , σ zr } are non-zero force-stress components; {m rθ , m θ r , m zθ , m θ z } are non-zero couplestress components; {ε rr , ε θθ , ε zz , ε rz , ε zr } are non-zero components of an infinitesimal strain tensor, {u r , u z } are non-zero components of the displacement vector; � θ is a non-zero component of the rotation tensor, {κ rθ , κ θr , κ zθ } are non-zero components of the curvature tensor; and µ are Lamé constants defined in the same fashion as that in the classical linear elasticity; and η and η ′ denote the material constants accounting for the presence of couple stresses. It is worth noting that η and η ′ are additional material parameters responsible for the length-scale effect (i.e., the presence of material microstructure) and, if these constants vanish, the couple stress theory will reduce identically to the classical linear elasticity.
Field equations for surface part. A material surface adhered to the top of the bulk is modeled by the surface elasticity theory proposed by Gurtin and Murdoch 26 , Gurtin and Murdoch 27 , and Gurtin et al. 28 . For an axisymmetric case, the non-zero surface displacements {u s r , u s z } , the non-zero surface strains {ε s rr , ε s θθ } , and the non-zero surface stresses {σ s rr , σ s θθ , σ s rz } are governed by where τ s denotes the residual surface tension; s , µ s are surface Lame's constants; and t s r , t s z are radial and vertical traction acting to the surface part by the bulk layer. Combining Eqs. (4)-(6) yields the equilibrium equations in terms of the surface displacements u s r , u s z as where κ s = 2µ s + s and the fact that the residual surface tension τ s is spatially independent has been utilized.
Boundary and continuity conditions. Since the surface part is perfectly bonded to the bulk layer, the surface displacements {u s r , u s z } and the tractions {t s r , t s z } can be related to the displacements and stress components of the bulk layer by www.nature.com/scientificreports/ By applying the continuity conditions Eqs. (9) and (10) together with surface Eqs. (7) and (8), it leads to a set of nonclassical boundary conditions on the top surface of the bulk layer: Since the surface part is considered infinitesimally thin and has no bending resistance, the applied couple traction m(r) on the top of the surface-bulk system is transmitted to the bulk layer directly and this yields an additional boundary condition: The boundary conditions at the bottom of the bulk layer can be readily expressed as Equations (11)-(16) form a complete set of boundary conditions for the bulk layer accounting for the surface effects.

Solution procedure
To obtain the closed-form solution of an elastic field within the bulk layer, a method of Hankel transform together with the representation of the displacement field is adopted. In particular, the vertical and radial displacements of the bulk layer undergoing the axisymmetric deformation admit the following representations 66,67 where α = ( + µ)/2( + 2µ) ; ℓ = √ η/µ represents the length scale of the bulk material; is an axisymmetric Laplacian operator; � = �(r, z) and � = �(r, z) are both solutions of the following equation: The closed-form general solution of Eq. (19) can be readily established by applying Hankel transform method 54,65,68 and the final results are given by where J m denotes Bessel function of the first kind of order m; ξ ∈ [0, ∞) is a transform parameter; ζ = 1 + ℓ 2 ξ 2 ; and C i (i = 1, 2, ..., 6) are unknown coefficients. The general solutions of the displacements {u r , u z } , the rotation � θ , the force stress components {σ rr , σ θθ , σ zz , σ rz , σ zr } , and the couple stress components {m rθ , m θr , m zθ , m θz } can be obtained upon substitution of Eqs. (20) and (21) into Eqs. (17), (18), (2), and (3). The explicit expressions for the complete elastic field within the bulk layer, in terms of the unknown coefficients C i (i = 1, 2, ..., 6) , are reported in Supplementary Appendix for the sake of brevity.
By enforcing the boundary conditions given by Eqs. (11)-(16) together with the general solutions for {u r , u z , σ zz , σ zr , m zθ } given in Supplementary Appendix, it yields the following system of linear algebraic equations for determining the unknown coefficients C i (i = 1, 2, ..., 6) : www.nature.com/scientificreports/ where C = { C 1 C 2 · · · C 6 } T and the coefficient matrix A(ξ ) and the vector F(ξ ) are given explicitly by The solution of the system (Eq. 22) for each ξ ∈ [0, ∞) can be obtained numerically via standard linear solvers. Once C i (i = 1, 2, ..., 6) are solved, the elastic field within the bulk layer can be obtained from supplementary Eqs. (A1)-(A12). To evaluate all involved improper integrals, an efficient quadrature rule similar to that employed by Rungamornrat et al. 54 and Lawongkerd et al. 65 is adopted.

Results and discussion
Computed results for certain cases are first compared with existing benchmark solutions to verify both the formulation and solution procedure. The influence of surface and couple stresses on the elastic field within a thin material layer under various surface loads is subsequently investigated. To clearly demonstrate the individual and simultaneous effects on the size-dependent characteristics, results for four different models (i.e., Model-1 with both surface and couple stress effects, Model-2 with only surface effect (i.e., ℓ → 0 ), Model-3 with only couple stress effect (i.e., τ s , κ s → 0 ), and Model-4 without surface and couple stress effects (i.e., τ s , κ s , ℓ → 0 )) are reported and compared. For convenience in simulations and presentation of results, following normalized coordinates and parameters r = r/� , z = z/� , a = a/� , h = h/� , τ s = τ s /2µ� , and l 0 = ℓ/� with � = κ s /2µ denoting the length scale of the material surface are introduced.
Verification. In the numerical study, material parameters reported by Miller and Shenoy 48 and Shenoy 49 are employed. In particular, Lamé constants of the bulk material are taken as = 58.17 × 10 9 N/m 2 and µ = 26.13 × 10 9 N/m 2 , whereas surface Lamé constants and the residual surface tension are taken as s = 6.8511 N/m , µ s = −0.376 N/m , and τ s = 1 N/m , respectively. Consider first an elastic half space subjected to a uniformly distributed normal traction p 0 over a circular region of radius a as illustrated in Fig. 2a. To simulate the half-space medium within the current setting, the thickness of the layer is taken to be sufficiently large in comparison with a and the ratio h/a = 1000 is considered in the analysis. Results for the force stress component σ zz and the couple stress component m θ r versus the ratio a/ℓ are compared with those reported by Lawongkerd et al. 65 in Fig. 3 for z/a = 0.25 , r/a = 0.5 , and l 0 = 1 . It is seen that the computed results are in excellent agreement with the benchmark solutions for all four models.
Another verification is carried out for an elastic layer under a uniformly distributed normal traction p 0 acting on a circular region of radius a shown in Fig. 2b for the load Case A. Results for this particular problem were reported by Rungamornrat et al. 54 for the classical case and the case with only surface stress effect. To simulate these two special cases, the parameters τ s , κ s , ℓ and ℓ are taken to be sufficiently small for each scenario. The computed surface displacements (i.e., z = 0 ) are reported in Fig. 4 for a = 10 and h = 10 and the stress components at the normalized depth z = 0.25 are shown in Fig. 5 for a = 1 and h = 10 . The good agreement between the two sets of results additionally confirms the validity of the proposed scheme and derived solutions.
Influence of surface and couple stresses. In this section, results from a parametric study are reported to demonstrate the role of both surface and couple stresses on the predicted response and size-dependent behavior of a substrate coated by a thin coating layer under surface loads. In particular, the load transferring characteristics from the coating surface to the substrate and the influence of the coating-layer thickness are of primary interest. To also explore the influence of applied loads and their distribution, a coated system subjected to four representative surface loads acting on a circular region of radius a shown in Fig. 2b,c (i.e., Case A for a uniformly distributed normal traction p(r) = p 0 , Case B for the Hertzian normal traction p(r) = p 0 1 − (r/a) 2 , Case C for a linearly distributed radial shear traction q(r) = q 0 r/a , and Case D for a quadratically distributed radial shear traction q(r) = q 0 (r/a) 2 ) are considered. In simulations, the following material parameters E = 76 GPa , ν = 0.3 , κ s = 1.22 N/m , and τ s = 0.89 N/m 48,49 are utilized unless stated otherwise. Note in addition that only the case of comparable surface and couple stress effects is investigated and, to simulate such scenario, the two (23)  www.nature.com/scientificreports/ material length scales ℓ, � are taken as l 0 = 1 . The full discussion on the two size effects for a wide range of the ratio ℓ/� can be found in the work of Le et al. 63 and Lawongkerd et al. 65 .
To demonstrate the intensity of load transferring to the coated substrate, the vertical stress σ zz for the load Case A and load Case B and the shear stress σ zr for the load Case C and load Case D are reported in Fig. 6 for z/a = 1 , h/a = 1 , and a/ℓ ∈ {0.01, 1, 100} . Three values of the ratio a/ℓ are considered to represent cases when the size of a loading region (representing the external length scale) is much less than, comparable to, and much larger than the two material length scales. Note that both vertical and shear stresses are normalized by the maximum intensity of the applied surface loads to clearly observe the role of the coating layer in the reduction of the transferring stresses to the substrate. For the first two load cases (i.e., load Case A and load Case B), the normalized vertical stress attains the maximum magnitude at the center of the loading region and monotonically decays to zero as r/a increases for all models and values of a/ℓ (see Fig. 6a-c). As the size of the loading region becomes comparable to both bulk and surface material length scales, transferring vertical stresses to the substrate are clearly different for all four models (see Fig. 6b). Such finding confirms the important role of both surface and couple stresses when a fall within the range of ℓ, � . Clearly, the Model-2 and Model-3 cannot be used as the replacement of the Model-1. In addition, the presence of surface and couple stress effects clearly reduce the maximum transferring stress to the substrate in comparison with the classical case; in particular, the Model-1 yields the least value of the maximum transferring stress. When a is much less than ℓ, � (see Fig. 6a), the Model-2 and Model-3 still predict responses differently from the classical case, but the effect of surface stresses is more pronounced than that of the couple stresses. The transferring vertical stresses obtained from the Model-1 and Model-2 are comparable but very different from those from the Model-3 and Model-4. These results suggest that the Model-2 can be used in lieu of the Model-1 to simplify the calculations when a ≪ ℓ ∼ � . When a is much larger than ℓ, � (see Fig. 6c), the surface and couple stresses play an insignificant role in the predicted response; in particular, results from the Model-1, Model-2, and Model-3 are almost identical to the classical solutions. For this range of external and material length scales, the Model-4 is considered sufficient for simulating the response of interest. It is worth noting that changing the distribution of applied normal tractions does not alter the response characteristics except for the difference in magnitude resulting from the difference in the traction resultant. For the load Case C and load Case D, the magnitude of the normalized shear stress σ zr /q 0 transferring to the substrate increases from zero at the center of the loading region to its maximum at r/a ∈ [0.5, 1] and then decays asymptotically to zero as r increases (see Fig. 6d-f). It is worth pointing out that for these loading conditions, the role of the surface stresses on the maximum transferring shear stress is opposite to that of the couple stresses. Specifically, the surface stresses (the Model-2) tend to lower the maximum transferring shear stress from the classical case while the couple stresses clearly boost such maximum and also switch the direction of the shear stress. By comparing results for three different values of a/ℓ and two different distributions of the applied shear loads, a similar conclusion to the load Case A and load Case B can be drawn. In particular, as the size of the loading region reduces to be comparable to (or much less than) the two length scales ℓ, � , the Model-1 (or the Model-1 and Model-2) must be used to capture the size effects. Note also that the presence of both surface and couple stresses can either reduce (as a ≪ ℓ ∼ � ) or boost (as a ∼ ℓ ∼ � ) the maximum transferring shear stress to the substrate from the classical case.
Through-the-thickness profiles of the vertical stress σ zz at r/a = 0 for the load Case A and load Case B and the shear stress σ zr at r/a = 0.7 for the load Case C and load Case D are also reported in Fig. 7 for h/a = 1 and a/ℓ ∈ {0.01, 1, 100} . The specific values of r/a used to collect those results are associated with the location where the transferring stress to the substrate attains (for the load Case A and load Case B) or approximately attains www.nature.com/scientificreports/ (for the load Case C and load Case D) its maximum. It is evident from Fig. 7a-c that the vertical stress predicted by the Model-1, Model-2 and Model-3 decreases faster than those in the classical case as the depth z increases. Nevertheless, for load Case C and load Case D (see Fig. 7d-f), the Model-2 tends to boost the decay of the shear stress across the thickness of the coating layer from the classical case, but the Model-3 seems to lower such decay. The Model-1 accounting for both effects can either lower (see Fig. 7e) or boost (see Fig. 7d) the decay depending on the ratio a/ℓ . For all load cases considered, the through-the-thickness profiles of both vertical and radial shear stresses are strongly dependent on both surface and couple stress effects when a is comparable to (see Fig. 7b,e) or much less than (see Fig. 7a,d) ℓ, � and for the latter case, the surface effect is found to be more pronounced. Note in addition that changing the distribution of applied surface loads does not alter the trend of the predicted response.
To further illustrate the influence of the coating-layer thickness on the reduction of the transferring stress on the substrate when both surface and couple stress effects are present, the transferring vertical and shear stresses for the case of applied normal and radial shear stresses, obtained from the Model-1, are reported as a function of the normalized thickness h/a in Fig. 8 for a/ℓ ∈ {0.01, 1, 100} . Since the role of surface and couple stress effects for different surface load distributions are similar, the load Case A and load Case C are chosen as representative load cases for applied normal and radial shear tractions, respectively. For the load Case A, the transferring vertical stress σ zz is reported at r/a = 0 where it attains the maximum (see Fig. 6a-c). For this loading case, the increase in the thickness of the coating layer can significantly lower the maximum transferring vertical stress to the substrate for both the classical model and Model-1. However, the presence of both couple and surface stresses www.nature.com/scientificreports/ renders such reduction more pronounced when the size of the loading region either falls within the range of or is much smaller than the material length scales ℓ, � (see Fig. 8a). For the latter case (as a ≪ ℓ ∼ � ), the surface effect is the key responsible for such substantial reduction from the classical case. For the load Case C, it is chosen, for convenience, to report the transferring radial shear stress σ zr at r/a = 0.7 since the exact location of the maximum r/a varies from 0.5 to 1 (see Fig. 6d-f). It is seen from this set of results that while the transferring shear stress to the substrate decreases monotonically and asymptotically to zero as the coating-layer thickness www.nature.com/scientificreports/ increases when both surface and couple stresses are taken into account, but the presence of such effects can either enhance (as a ∼ ℓ ∼ � ) or reduce ( a ≪ ℓ ∼ � ) the transferring shear stress from the classical case. In addition, the switch of the direction of the transferring shear stress from that of the applied shear traction for a sufficiently large h/a , as observed for the classical case, disappears when both surface and couple stress effects are significant. Finally, the size dependent characteristics of the predicted transferring stresses to the substrate are also investigated. To clearly illustrate such behavior, the maximum transferring vertical stress for the load Case A www.nature.com/scientificreports/ and the transferring shear stress at r/a = 0.7 for the load Case C are reported as a function of a/ℓ in Fig. 9 for h/a ∈ {0.5, 1, 2} . It is seen that for any given aspect ratio h/a , the normalized transferring stresses to the substrate obtained from the Model-1 are strongly size dependent or, equivalently, depend on the length scale ratio a/ℓ . As a decreases to be comparable to or less than ℓ , the maximum transferring vertical stress for the load Case A drops quite rapidly and monotonically from the value predicted by the classical model. The different behavior is observed for the load Case C. The variation of the transferring shear stress over a wide range of the ratio a/ℓ is not monotone; in particular, for a comparable to or larger than ℓ, � , the predicted transferring shear stress from the Model-1 is higher than the classical solution while the reverse trend can be concluded when a is much less than ℓ, � . For both loading cases, the size dependency decays to be insignificant as the size of the loading region a is much larger than the material length scales ℓ, � , and the Model-4 is therefore sufficient for the simulations.

Conclusion
The analytical solution of an elastic field of a thin material layer coating a rigid substrate and excited by axisymmetrically distributed surface loads have been derived. Such established results are considered novel in that both surface free energy and material microstructures, which are recognized to be responsible for the size effects in  www.nature.com/scientificreports/ small-scale objects, are taken into account simultaneously to model a medium of finite thickness. This allows the direct application to simulate the mechanical response of components coated by a very thin material layer.
A continuum-based model integrating the couple stress elasticity theory for tackling the inherent microstructural effect and Gurtin-Murdoch surface elasticity theory for capturing the surface effect has been formulated and solved by an analytical scheme based on Hankel transform and the displacement representation. Obtained results are explicit in an integral form, highly accurate as useful benchmark solutions, and an essential basis for the development of solution schemes to tackle surface contact problems.
Results from an extensive numerical study have revealed that the surface and couple stresses significantly affect both the characteristics and maximum value of load transferring to the coated substrate in comparison with the classical case when the size of the loading region is comparable to or much smaller than the material length scales. For a coated system under normal tractions, the presence of surface and couple stresses can significantly boost the reduction of the transferring vertical stress to the substrate especially when the size of the loading region is much less than the length scale of the bulk and surface materials. A different trend has been observed in the case of applied shear loads. The transferring shear stress to the substrate predicted by the model integrating both surface and couple stress effects can be either lower or higher than the classical solution depending on the ratio between the size of the loading region and the material length scales. This results directly from that the surface effect lowers the transferring shear stress but the couple stress effect causes the reverse trend. In addition, as the size of the loading region become much smaller than the two material length scales, the surface effect is much more pronounced than the couple stress effect.

Data availability
The datasets generated and/or analyzed during the current study are not publicly available due to that the data also forms part of an ongoing study, but are available from the corresponding author on reasonable request.